R Notebook

library(tidyverse)
library(lubridate)
library(tsibble)
library(tsibbledata)
library(fable)
library(nycflights13)

1 Question 1 Load in the nyc_bikes data from the tsibbledata package. Have an initial look at it to see what you’re working with. Create three new columns: one that stores only the year the bike was used, one that stores only the month the bike was used, and one that stores the date. Use the data stored in start_time to create these new columns.

nyc_bikes

codeclan solution

nyc_bikes_year <- nyc_bikes %>%
  select(-gender, -type, -birth_year) %>%
  mutate(year_of_use = year(start_time),
         month_of_use = month(start_time, label = TRUE), 
         date_of_use = date(start_time)
  )

head(nyc_bikes_year,3)

NA

2 Question 2 Summarise the number of bike hire counts by month. Make a plot of this data. *Hint: remember that to group time series (tsibble) data, you need to use index_by instead of group_by before your summarise function. What does this plot tell you about the time series? Do you think this downsampled data would be adequate to build a forecast with?

nyc_count <- nyc_bikes_time %>%
  group_by_key(bike_id) %>%
  index_by(month) %>%
  summarise(bike_hire_counts = n())

ggplot(data = nyc_count, aes(x = month, y = bike_hire_counts, group = bike_id)) +
  geom_line(aes(color=bike_id))

Codeclan Solution

nyc_bikes_month_summary <- nyc_bikes_year %>% 
  index_by(month_of_use) %>%
  summarise(bike_hire_counts = n())

ggplot(nyc_bikes_month_summary) + 
  aes(x = month_of_use, y = bike_hire_counts) +
  geom_point() + 
  geom_line(group = 1) +              
  ggtitle("Bike hire in NYC") + 
  xlab("Month") + ylab ("total count")

3 Question 3 Now Summarise the number of bike hire counts by date. Make a plot of this new aggregated data. What does this plot tell you about the time series? Would this data be preferrable for time series forecasting compared to the monthly data?

nyc_bikes_date <- nyc_bikes_time %>% 
  index_by(date) %>%
  summarise(bike_hire_counts = n())

nyc_bikes_date 

ggplot(data = nyc_bikes_date, aes(x = date, y = bike_hire_counts)) +
  geom_line()

NA
NA

CODECLAN SOLUTION

# count by date 
nyc_bikes_date_summary <- nyc_bikes_year %>% 
  index_by(date_of_use) %>%
  summarise(bike_hire_counts = n())

nyc_bikes_date_summary

ggplot(nyc_bikes_date_summary) + 
  aes(x = date_of_use, y = bike_hire_counts) +
  geom_point() + 
  geom_line(group = 1) +              
  ggtitle("Bike hire in NYC") + 
  xlab("Date") + ylab ("total count")

Let’s go with our date one. The reasons for this is as follows:

1. Forecasts work better with more data. The monthly data didn’t have that many data points.
2. There will be daily variation, and monthly variation that will be more accurately captured by using date based data.

4 Question 4

Let’s begin to build a model. For this, we will use the downsampled by date dataset we created above in question 3. If you haven’t managed to get the code to do this, you can find it below.

count <- nyc_bikes_date %>% 
  select(date, bike_hire_counts)

count %>%
  autoplot()

Plot variable not specified, automatically selected `.vars = bike_hire_counts`

fit_bikes <- nyc_bikes_date %>%
    fill_gaps() %>%
    model(
    snaive = SNAIVE(count),
    mean_model = MEAN(count),
    arima = ARIMA(count)
  )

argument is not numeric or logical: returning NA1 error encountered for snaive
[1] invalid time series parameters specified
1 error encountered for mean_model
[1] (list) object cannot be coerced to type 'double'
1 error encountered for arima
[1] 'x' must be a vector of an atomic type

# fill the gaps : here i'll use median as it's a decent summary of how bike rentals work (not overly affected by seasonal changes)
nyc_bikes_filled <- nyc_bikes_date %>%
  fill_gaps(bike_hire_counts = as.integer(median(bike_hire_counts)))

fit_bikes <- nyc_bikes_filled %>%
    model(
    snaive = SNAIVE(count),
    mean_model = MEAN(count),
    arima = ARIMA(count)
  )

argument is not numeric or logical: returning NA1 error encountered for snaive
[1] invalid time series parameters specified
1 error encountered for mean_model
[1] (list) object cannot be coerced to type 'double'
1 error encountered for arima
[1] 'x' must be a vector of an atomic type

CODECLAN SOLUTION

# count by date 
nyc_bikes_date_summary <- nyc_bikes_year %>% 
  index_by(date_of_use) %>%
  summarise(bike_hire_counts = n())

nyc_bikes_date_summary

# fill the gaps : here i'll use median as it's a decent summary of how bike rentals work (not overly affected by seasonal changes)
nyc_bikes_filled <- nyc_bikes_date_summary %>%
  fill_gaps(bike_hire_counts = as.integer(median(bike_hire_counts)))

# build a model
bikes_fit <- nyc_bikes_filled %>%
  model(
    naive = NAIVE(bike_hire_counts), 
    mean_model = MEAN(bike_hire_counts), 
    snaive = SNAIVE(bike_hire_counts)
  )

5 Question 5 Now we have our model fit, build a forecast to predict bike use over the next four months. Plot your models alongside your data. Hint: forecast parameter would be roughly 120 (30 days x 4 months)

# build a forecast
bikes_forecast <- bikes_fit %>%
  fabletools::forecast(h = 30*4)
bikes_forecast

bikes_forecast %>%
  autoplot(nyc_bikes_filled, level = NULL) 

# create a shorter interval to plot: years after 1980
bikes_shorter <- nyc_bikes_filled %>%
  mutate(month = month(date_of_use)) %>%
  filter(month >= 8)

# check one model at a time
bikes_forecast %>%
  filter(.model == "snaive") %>%
  autoplot(bikes_shorter, level = NULL) 

6 Question 6 Test your model accuracy : choose a training data set from your main dataset, build a forecast on the training set, and then plot the training set forecast against the real data. Calculate model accuracy.

## test model accuracy 
# Set training data from Jan to Nov
train <- nyc_bikes_filled %>%
  mutate(month = month(date_of_use)) %>%
  filter(month <= 10)

# run the model on the training set 
bikes_fit_test <- train %>%
  model(
    naive = NAIVE(bike_hire_counts), 
    mean_model = MEAN(bike_hire_counts), 
    snaive = SNAIVE(bike_hire_counts))

bikes_fit_test

# forecast from the training set
bikes_forecast_test <- bikes_fit_test %>%
  fabletools::forecast(h = 12*6)

# Plot forecasts against actual values
bikes_forecast_test %>%
  autoplot(train, level = NULL) + 
  autolayer(nyc_bikes_filled, color = "black")

Plot variable not specified, automatically selected `.vars = bike_hire_counts`

# note accuracy, arrange by best
accuracy_model <- fabletools::accuracy(bikes_forecast_test, nyc_bikes_filled)

The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
11 observations are missing between 2019-01-01 and 2019-01-11

accuracy_model %>% 
  select(-.type) %>%
  arrange(RMSE)

NA

7 Question 7 Look at your forecast plots and accuracy values. Describe your results. Are your models a good fit for the data? If not, why not? What would you suggest doing with the data if you were expected to present these back to a client? For example, would you ask for more data? Would you test a different model?

SUMMARISE YOUR FINDINGS

Forecast isn’t great. This could be for several reasons:

1. not enough data
2. something is affecting it
3. wrong model choice
4. I’d collect more data if possible - we only have a year’s worth so far.

8 Question 8 Make a simple ggplot (geom_point) which plots the start longitude and latitudes of each bike. Create a separate facet for each bike_id. Colour the dots in by month of use. What does this tell you about what month each bike was used most in?

Do the same for the end longitude and latitudes.

# make the month 
nyc_bikes_spatial <-nyc_bikes %>%
  select(bike_id, start_time, start_lat, start_long, end_lat, end_long) %>%
  mutate(month = month(start_time, label = TRUE))

# make a faceted plot : start points
ggplot(nyc_bikes_spatial) + 
  geom_point(aes(x=start_lat, y = start_long, colour = month)) +
  facet_wrap(~bike_id) +
  ggtitle("NYC bikes: start points") +
  xlab("Latitude") + ylab("Longitude")

# make a faceted plot : end points
ggplot(nyc_bikes_spatial) + 
  geom_point(aes(x=end_lat, y = end_long, colour = month)) +
  facet_wrap(~bike_id) +
  ggtitle("NYC bikes : end points") +
  xlab("Latitude") + ylab("Longitude")

9 Question 9 Create an interactive leaflet plot which plots the start points of the city bikes. Ensure it has at least markers to denote start points (taken from the nyc_bikes_spatial data). Feel free to add any additional features you wish.

# make a leaflet plot
library(leaflet)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

# create leaflet map
leaflet(nyc_bikes_spatial) %>%
  addTiles() %>%
  addCircleMarkers(lng = ~start_long,
                   lat = ~start_lat) %>%
    setView(-74.00, 40.71, zoom = 12) %>%
  addProviderTiles("CartoDB.Positron")

NA

10 Extension This is VERY much an extension task, and will take a bit of time (as well as googling) if you decide to attempt it.

Adapt your leaflet plot which adds lines showing the start and end point of each bike

install.packages("geosphere")

trying URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.6/geosphere_1.5-10.tgz'
Content type 'application/x-gzip' length 1003838 bytes (980 KB)
==================================================
downloaded 980 KB

The downloaded binary packages are in
    /var/folders/58/m5fvfpw93rz5rtzg6nc5_lbr0000gn/T//Rtmpl1CVEN/downloaded_packages

# load the geosphere package to add intermediate points (lines)
library(geosphere)

# subset the data to make it easier
bikes <- nyc_bikes %>%
  select(bike_id, start_time, start_lat, start_long, end_lat, end_long)

# get start points (longitude has to go first)
start <- bikes %>%
  select(start_long, start_lat)

end <- bikes %>% 
  select(end_long, end_lat)

# convert to matrix  
start_matrix <- as.matrix(start[ , c(1,2)])
end_matrix <- as.matrix(end[ , c(1,2)])

# plot using leaflet
gcIntermediate(start_matrix, end_matrix,  
           n=100, 
           addStartEnd=TRUE,
           sp=TRUE) %>% 
leaflet() %>% 
addTiles() %>% 
addPolylines() 

NA